
DEPARTMENT OF MATHEMATICAL SCIENCES 
COLLEGE OF SCIENCES 
OLD DOMINION UNIVERSITY 
NORFOLK, VIRGINIA 23529 




FINITE— VOLUME APPLICATION OF HIGH ORDER ENO SCHEMES 
TO TWO-DIMENSIONAL BOUNDARY-VALUE PROBLEMS 


By 

J. Mark Dorrepaal, Principal Investigator 
and 

Jay Casper, Graduate Research Assistant 
Final Report 

For the period ended August 15, 1990 



f 


Prepared for . . . 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, Virginia 23665 


Under 

Master Contract Agreement NAS1-18584 
Task Authorization No. 56 
James D. Keller, Technical Monitor 
FLDMD-Theoretical Aerodynamics Branch 


N ASA— CP-1 86917) PINT T e- VUL U M t APPLICATION 
)F HIGH ORDER LNG SCHEMES TO 
IULT I— DIMENSIONAL BOUNDARY- VALUE PROtfLLMS 
i na 1 Report, period ending 15 Aun. 1990 
[Old Dominion Univ.) J9 r> CSCL 12A G3/oA 


N90-28387 


Unc 1 as 
0302492 


August 1990 



DEPARTMENT OF MATHEMATICAL SCIENCES 
COLLEGE OF SCIENCES 
OLD DOMINION UNIVERSITY 
NORFOLK, VIRGINIA 23529 


FINITE— VOLUME APPLICATION OF HIGH ORDER ENO SCHEMES 
TO TWO-DIMENSIONAL BOUNDARY -VALUE PROBLEMS 


By 

J. Mark Dorrepaal, Principal Investigator 


and 

Jay Casper, Graduate Research Assistant 


Final Report 

For the period ended August 15, 1990 


Prepared for . , . 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, Virginia 23665 


Under 

Master Contract Agreement NAS1-18584 
Task Authorization No. 56 
James D. Keller, Technical Monitor 
FLDMD-Theoretical Aerodynamics Branch 


Submitted by the 

Old Dominion University Research Foundation 
P.O. Box 6369 

Norfolk, Virginia 23508-0369 


August 1990 


ABSTRACT 


Finite- Volume Application of High Order ENO Schemes 
to Multi-Dimensional Boundary-Value Problems 


Jay Casper 

Department of Mathematics 
Old Dominion University 
Norfolk, Virginia 


In recent years, a class of numerical schemes for solving hyperbolic partial differential 
equations has been developed which generalizes the first-order method of Godunov to arbi- 
trary order of accuracy. High-order accuracy is obtained, wherever the solution is smooth, 
by an essentially non-oscillatory (ENO) piecewise polynomial reconstruction procedure, 
which yields high-order pointwise information from the cell averages of the solution at 
a given point in time. When applied to piecewise smooth initial data, this reconstruc- 
tion enables a flux computation that provides a time update of the solution which is of 
high-order accuracy, wherever the function is smooth, and avoids a Gibbs phenomenon at 

discontinuities. „ , . . , . . , . , , u . 

The application of these schemes to areas of scientific and industrial interest obvi- 
ously requires compressible flow solutions in more than one spatial dimension. The multi— 
dimensional ENO schemes proposed thus far in the literature are of finite-difference type, 
combining an adaptive high-order ENO interpolation procedure with a split-flux approach. 
Due to the design of these schemes, it is not clear, at present, how to apply them at bound- 
aries, t.g. solid walls. Moreover, as with any finite-difference scheme in which stencils are 
not pre-determined, curvature and discontinuity of the computational grid can pose fun- 
damental problems. , . , ... ,. . i 

In this paper, we consider the finite-volume approach in developing multi-dimensional, 
high-order accurate ENO schemes. In particular, a two-dimensional extension is proposed 
for the Euler equations of gas dynamics. This requires a spatial reconstruction operator 
that attains formal high order of accuracy in two dimensions by taking account of cross 
gradients. Given a set of cell averages in two spatial variables, polynomial interpolation of 
a two-dimensional primitive function is employed in order to extract high-order pointwise 
values on cell interfaces. These points are appropriately chosen so that correspondingly 
high-order flux integrals are obtained through each interface by quadrature, at each point 
having calculated a flux contribution in an upwind fashion. The solution-m-the-small of 
Riemann’s TVP that is required for this pointwise flux computation is achieved using Roe s 
approximate Riemann solver. Issues to be considered in this two-dimensional extension 
include the implementation of boundary conditions and application to general curvilinear 

coordinates. , . . 

Results of numerical experiments are presented for qualitative and quantitative exami- 
nation. These results contain the first successful application of ENO schemes to boundary- 
value problems with solid walls. 


Finite- Volume Application of High Order ENO Schemes 
to Two-Dimensional Boundary- Value problems 


Jay Casper 


Department of Mathematics 
Old Dominion University 
Norfolk, Virginia 


l. Introduction 

We wish to design high-order accurate essentially non-oscillatory (ENO) schemes for 
the numerical approximation of weak solutions of hyperbolic systems of conservation laws 

u t + /(u)* + y(u)„ = 0 , ( la ) 


subject to given initial conditions 

u(i,y,0) = u°(i, y) . 


(lb) 


The function u = (u 1 ^ 2 , . . . ,u m ) T is a state vector and the fluxes /(u) and y(u) are 
vector-valued differentiable functions of m components. We assume that the system (la) 
is hyperbolic in the sense that the m x m Jacobian matrices 



»<“> = Tu 


are such that any linear combination of A and B has m real eigenvalues {A (tt)} and a 
complete set of m right eigenvectors {r k (u)> and left eigenvectors {/*(«)}, which we assume 
to satisfy the orthonormal relation 1’ • r ] = <5,; . 
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We assume that the initial-value problem (l) is well-posed in the sense that the solution 
u depends continuously on the initial data, and that this solution is generically piecewise 
smooth, with at most a finite number of discontinuities. Appealing to well-known theoret- 
ical results, we know that a weak solution of the IVP (1) must satisfy (la) in the sense of 
distribution theory. Furthermore, this is equivalent to requiring that u obey the integral 
form of (la), where the limits of integration can reflect any smoothly bounded domain in 
the x— y plane and any time interval (t*,t 2 ) ■ Seeking such a solution, let 


x,-_ 1/2 < x < x i+ 1 / 2 , yj -i /2 < y < yj+ 1/2 , -» < *\i < 00 > 


denote a rectangular partition of the x — y plane, with x,- and y, denoting the centroids 
of their respective intervals. With a semi-discrete formulation in mind, we note that, for 
every every rectangle (xj_ x / 2 , £1+1/2) x (y*-i/2> Vj+i/i) > a weak solution of (la) must satisfy 


Q 1 

— a„(t) = fi+l/2 ,* (0 - fi— 1/2 ,2 (0 + &.J+ 1/2(0 - 

at a,ij 

where a,-,- is the area of the rectangle and 


(2a) 


«<,( 0 = — f + /a [ V + 1/1 u ( x ,y,t) dxdy (2b) 

&ij * **—1/2 ^ Vy-1/2 

is the cell average of u over the control volume at time t. The fluxes f and g are given by 


A 

/ i+1/2 ,) (0 = 

f Vi+l/i /(u(x <+ 1 / 2 , y, t))dy ; 

(2c) 


'V/-J /9 


&, *+1/2(0 = 

f 1+1/3 g{u{x, yj+ 1/2 , t))dx . 

(2d) 


We now treat (2a) as a system of ordinary differential equations for the purpose of 
time discretization, using a “method-of-lines” approach. Along any t = constant line, the 
right-hand side of (2a) is strictly a spatial operation in the unknown u, and we rewrite 
this equation, for fixed t, in the abstract operator-product form 

= ( £ «(() )„ , ( 3 ) 

thus effectively “separating” the spatial and temporal operations for computing solutions 
of (1). 

In [8,9], Runge-Kutta methods are presented for the time discretization of ordinary 
differential equations which are of high-order accuracy and TVD, in the sense that the 
temporal operator itself does not increase the total variation of the solution. We will employ 
these Runge-Kutta methods in order to achieve our desired accuracy in time. The use of 
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the formulation (2a), rather than the finite-difference approach of [8,9], makes for a more 
immediate application of its numerical approximation to the solutions of boundary-value 
problems involving solid walls or non-trivial geometries. 


2. Two-Dimensional ENO Schemes 

Given {&”} , cell averages of a piecewise smooth solution u(x, y, t”) of (1), we desire a 
high— order accurate numerical solution operator which will update these averages to 
time t n+l = t n + At. Specifically, we require that E h be r th -order accurate in the sense of 
local truncation error, i.e. when applied to the exact solution, E^ satisfies 

E h 6" - fl n+1 = 0(h r+l ) (4) 

in a pointwise sense, wherever u is sufficiently smooth, with Ax , Ay , and At all assumed to 
be 0(h) . Furthermore, we desire that our numerical update scheme avoid the development 
of spurious 0(1) oscillations near discontinuities in u. In order to achieve this property we 
require our operator to be essentially non-oscillatory (ENO). Such an operator satisfies, 
in the one-dimensional scalar case, 

TV{ E h a) = TV (u) + 0(h 1+p ) , (5) 

for some p > 0 , where TV represents total variation in x. For details concerning the initial 
development of ENO schemes, the reader is referred to [5] and the references therein. 

Employing the formulation (3) for the numerical update of the {& } , we discretize the 
temporal operation by using a Runge-Kutta method, which we write in the form 

fi. ! ? = El**- Sir 1 + Am A< (£«'">)« ] , 1 = 1,2 P , («) 

m=0 



6i;> = a? 1 


The order of accuracy achieved by this time discretization, as well as its TVD property, 
is determined by the values of the integer p and the coefficients a and (3 . ( See [8,9] for 
details. ) 

Now, if we assume that the scheme (6) achieves our desired r th — order accuracy in 
time, then clearly, this scheme satisfies (4) if we can evaluate ( £ u(t) ),y , the exact spatial 
operation on the right-hand side of (2a). However, the calculation of the fluxes (2c— d) 
needed for this evaluation requires that we know the solution u(x,y,t), pointwise, at a 
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given time i, whereas the information we have at any time t is that of the cell averages 
(2b). And clearly, since 

MO = «*(*«. Vi. 0 + °( hi ) . 

wherever u is smooth, there is an inherent limit on the order of accuracy if we use the 
cell averages themselves in the flux calculation. Therefore, we replace the operator Z with 
a discrete spatial operator L which acts upon the (u(i)} and approximates the pointwise 
operation of Z to high order. To this end, we see that if 

L fl(t) = Z u{t) + 0{h r ) , (7a) 


wherever u is smooth, then, upon replacement of Z in (5), the local truncation error of 
our fully discrete scheme will satisfy (4), as required. 

We define L explicitly by 




_ 1 _ 

a <i l 


7i+ 1/2, i(0 


- 7i-l/2,;(*) + ft\i+l/*(0 




(7b) 


where J and <? are to be designed as high-order numerical approximations to f and g in 
(2c-d). We require that the numerical flux functions J , g 

Ji+i/tA*) = (0 «*+*.,•+-(*)) * 1 <k,m<r, 

(7c) 

5t,j+l/2(0 = 5(®»-*+l,i-n+l(0i* • * 

be Lipschitz continuous functions of their arguments and consistent with the true fluxes 
7(u) , g(u) in the sense 

J(w , ...,«;) = f{w), g(ti;,...,u>) = g{tv) . (7d) 

A . 

The first and most important step in the high-order approximations of / , g is the 
method by which we obtain high-order accurate pointwise information of the solution 
u(z, y,t) from the given set of cell averages (6(t)} . Let R* be a spatial operator which 
reconstructs this set of cell averages and yields a two-dimensional, piecewise polynomial 
R 2 {x, y ; fi(t)) of degree r-1 which approximates u(x, y, t) with a truncation error of 0(h r ) , 
wherever u is sufficiently smooth. We write this relationship in the form 

R 2 (x,y;U{t)) = u(x,y,t) +e(x,y)h r . (8) 

In general, the “reconstruction” operator R 2 which we have developed is a “natural” two- 
dimensional extension of the one-dimensional “reconstruction by primitive” presented in 
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[5], The reconstruction procedure involves polynomial interpolation in combination with 
an adaptive stencil algorithm which, in tandem, achieve high-order accuracy in smooth 
regions while avoiding 0(1) oscillations near steep gradients. This operator is detailed in 
Section 3. 

In order to include the more general case, where /(u) cannot be integrated in closed 
form in (2c), we will approximate this integration by Gaussian quadrature. In order 
to express the error made by this approximation, let q(x) be a function whose 

integration on [a, 6] we approximate by the “classical” Gaussian quadrature, i.e. relative 
to the unit weight function on the interval [—1,1] . It can be shown (e.g. [10]) that the 
error made by this approximation with a K-point quadrature is given by 



6 — a 
2 


£ 


Jb=l 


c» 9(*») 


(MO* 



dx , 


for some £ in (a, b ) , with P K being the polynomial of degree K in the orthogonal basis 
that spans the space of polynomials of degree not exceeding K . This quadrature is exact 
when y(x) is a polynomial of degree less than or equal to 2 K — 1 . The roots of Pk{ x ) — 
(x — Zi)(x— £ 2 ) ... (i-ijr) are real and distinct, making it clear that the above truncation 
error is 0{h 2K+1 ) . Relating this error to the (r-l)-st degree polynomial reconstruction 
(8), we see that for r < 2K , this truncation error is at worst 0(h r+1 ) when r-1 is odd, 
and 0{h T+ 2 ) when r — 1 is even. Therefore, using the “larger” error, for fixed x and t, and 
sufficiently smooth /, the approximation of the flux integral (2c) by Gaussian quadrature 

satisfies 


%+i/tA*) = r* /(«(**+!/> ’ 0)<*y 


= M £ c, /( «(*,+,/, ,«,!)) + *(*+./: . l) '‘ r+1 , (9a) 

^ Jk=l 


for some rj in (j/y_i/ 2 , t/y+ 1 / 2 ) • 

Let v*(x, y, t) denote the piecewise polynomial approximation to u which is determined 
by the reconstruction operator R 2 in (8) and therefore satisfies 

Vh(x,y,t) = u(x,y,t) +e(x,y)h r , (9b) 

for fixed t, wherever u is sufficiently smooth. Since / is assumed differentiable in u, it is 
therefore Lipschitz continuous in u, and thus, for fixed t, 

/(u(x, y,t)) = f{v h {x,y,t)) +d{x,y)h r , (9c) 
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where d(x,y) = 0(e(x,y)) . Finally, we substitute (9c) into the quadrature in (9a), and 
we see that 

= “k 53 Cfc ( / ( V h (x i+ i /2 , y t , 0 ) + d (*i+ 1/2 1 y*) h ' 1 ( 9d ) 

Z *=1 

+ s(x,+ 1/2 , r?) h r+1 . 

Therefore, if we define our “abstract” numerical flux 7«+i/2,j(0 i* 1 (7b) by 

W.,(<) = =r E «*/( .».<)). (“») 

Z *=1 

then the error made by the approximate flux difference 7<+i/a,i(0 ~ 7*-i/a,j(0 ^e 
definition (7b) is given by 

/i+l/2,i(0 — /»'— l/2,i(0 = 7t+l/2,j'(0 — 7 <—1/2 .j (0 

+ X) c* [ d(x i+ i / 2 , y*) - d(x,_ 1 / 2 , yt) ] h r 

1 jk=i 

+ [ s(x, + 1 /2 , T?) - «(x<_i/2 , »?) ] h T+l . (10b) 

Clearly, if d and s are Lipschitz continuous on [x,_i/ 2 ) ^i+ 1 / 2 ] f° r each y , then the error 
relation in 10(b) satisfies 

7i+l/2,j(0 /• — l/2,j(0 = 7i+l/2,i(0 7»— l/2,j (0 + 0(h r+ *) . (10°) 

Moreover, a symmetrical argument can be used to show 

9i,j+i/t{t) ~ £»,j-i/2(0 = 5 »,j'+i/ 2(0 — ?«,j-i/2(0 + 0(h + ), (10d) 

where 

?f,i+i/a(0 = X) c * &( . yi+i/* » 0 ) • ( 10e ) 

L k=l 

Noting that the area a, ; is 0(h 2 ) , we see that upon substitution of the numerical fluxes 
(10a) and (lOe) which satisfy the error relations (10c) and (10d) into (7b), we have thus 
designed the spatial operator L that satisfies (7a), and therefore, when substituted for £ 
in (6), yields a numerical solution operator which is r th - order accurate in the sense of 
local truncation error. We note here that the desired truncation errors given by (10c— d) 
are achieved only if the functions representing the errors due to the quadrature (9a) and 
the solution approximation (8) are Lipschitz continuous on [x<-i/ 2 , £*+ 1 / 2 ] x [s/j— 1 / 2 » yy+ 1 / 2 ] • 
We now wish to modify the “abstract” numerical fluxes (10a) and (10e) such that 
(4) still holds in regions where the solution is smooth and, in addition, these fluxes will 
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account for possible discontinuities in u . This modification is largely due to the nature 
of the reconstruction step, at which the pointwise behavior of u(i, y,t) is approximated 
in a piecewise polynomial fashion within each cell. (See Section 3.) In particular, the 
piecewise polynomial generated by the reconstruction operator can be discontinuous at 
cell interfaces, due to an adaptive interpolation procedure. The relative size of these local 

“jumps” is on the level of the interpolation error in smooth regions and is 0(1) near 

discontinuities in u. Therefore, in order to resolve these discontinuities, the flux integrands 
in (2c-d) are approximated by 

/(«(*.».*)) » /““(^(x-O.y; O(0)> R 2 {x + 0,y\ fi(i))] , (11a) 

y(u(x, y, 0) « </ Rm [ R 2 (x , y — 0; &(t) ) , R 2 ( x, y + 0 ; fi(t) ) ] , (lib) 


where / Rm [u!,u 2 ] denotes the flux across x = 0 of the solution to the Riemann problem 
whose initial states are u x and u 2 . We use the notation q[x + 0) , q(x - 0) to denote 
the limiting values of q at x from the right and left, respectively. When the solution u 
is sufficiently smooth, the “jumps” in the discontinuities in the approximate solution at 
cell interfaces will be 0(h r ) in which case our previous conclusions concerning high-order 
truncation error are not altered by the modification (11). 

In fully-discrete form, we write our scheme as the Runge-Kutta method 


&g> = E i = i,2 ,. -p, (i2») 


m=0 


«<"> = , 


i,W = 6" + ‘ , 


where 


and 


(£ *'“’)« = - - 


->(m) 7(m) , ^(m) _ «(m) 

f i+1/2 ,j “ * i—1/2 1/2 »*,J- 1/» 


(12b) 


Ck /Rmf Ri (z,+1/J ” 0 ’ Vk ; ® (m)) ’ R ^ Xi+l/2 + 0 ’ Vk ; 1 ’ (12c) 


a (m) 
*,i+ 1/* 




k=l 


Assuming the error functions d(x, y) and s(x,y) in (10b) to be globally Lipschitz con- 
tinuous, the numerical solution operator E h defined by (12) is formally r th -order accurate 
in the sense of local truncation error as given by (4). Furthermore, if these error functions 
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remain Lipschitz continuous for N time steps, where N = t/ At = 0(l/h ) , we assume 
the cumulative error to be 0(h r ) . Thus, at the end of such a computation, we have a set 
} , approximations to the cell averages of u at time t N which satisfy 

% N - U N = 0(h r ) . (13a) 

If we desire our high— order accurate output in pointwise form, we simply perform one final 
reconstruction which, by (8) and (13a), will yield 

R 2 (x,y-,V N ) = u(x,y,t N ) + 0(h r ) . (13b) 

In addition to the accuracy properties (13), we desire that if u should develop disconti- 
nuities, then the scheme (12) will avoid 0(1) spurious oscillations, and we will design the 
reconstruction operator R 2 to do so. 

In Section 7, we present results which support the claim of cumulative error represented 
by (13). We will also extend the scheme (12) to general curvilinear co-ordinates and 
perform numerical experiments which apply this scheme to boundary— value problems. 
Some of these results are shown in Section 7, and are representative of the first successful 
attempt to employ high-order ENO schemes to two-dimensional problems with solid walls. 


3. Reconstruction 


Before proceeding to a discussion of the implementation of the scheme (12) in solving 
boundary-value problems, we turn now to describe more fully the high-order spatial recon- 
struction operator R 2 which is crucial to the scheme’s accuracy. For the purpose of clarity, 
we discuss the finer points of this procedure within the framework a scalar function defined 
on a rectangular computational mesh. Furthermore, as there are many ways to approach 
the development of this high-order reconstruction, for simplicity, we choose to describe the 
implementation of R 2 as a composition of two applications of a one-dimensional operator 
R, where the latter is the “reconstruction-by-primitive” operator in [5]. 

In preparation, we find it necessary to briefly review the notion of essentially non- 
oscillatory interpolation. Let w(x) be a piecewise smooth function in x whose values are 
known to us only at the discrete points {x,}. We introduce Q n (x] w) , an n th -degree 
piecewise polynomial function of x that interpolates w at the points {x^} , i.t. 


Q n {x\xu) = q n i(x]w) , 


(14a) 

(14b) 
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Xi ^ X ^ X| + i j 


where q n ,i is the (unique) n th -degree polynomial that interpolates w(x) at n+1 successive 
points {xy} which include x, and x,+i. We are therefore free to chose the other n-1 points, 
and we do so subject to the condition that u;(x) be “smoothest” on the chosen stencil in 
some asymptotic sense. 

It can be argued that information relevant to the smoothness of tv cam be extracted 
from a table of divided differences of tv. Employing the standard notation, the k th divided 
difference of w can be defined recursively by 


w[xj ] = tv(xj) ; 


IV Xy , . . . , Xj+ic] 


U>[xy+i, . . . , Xy+fc] W Xy, . . . , 3<y+fc— 1, 

x j+k ~ x j 


In particular, it can be argued that the magnitude of a k th divided difference, provides 
an asymptotic measure of the smoothness of tv in (xy, Xy+*) , in the following sense. Suppose 
that tv is smooth in the interval (io,x*) but is discontinuous in (xi,x*+i). Then for h 
sufficiently small, we expect 


|u>[x 0 ,...,Xfc] | < |u;[x 1 ,...,x* +1 ]| , 

and hence these divided differences can serve as a tool to compare the relative smoothness 
of w in various stencils. Since we always assume any stencil we choose to be contiguous, 
we can assign a particular stencil by determining its left-most index, which we denote by 
j(t). We choose j(») with the following hierarchical algorithm. 

Let jifc(t) denote the left-most index of a chosen “smoothest” (fc+l)-point stencil which 
includes the interval (x t ,x,+i) . Denote this stencil 

{ x h(*)-> • • • > x it(0+ fc fc = l,2, ...,n (15a) 

Since any stencil must include {xj,x,+i} , our recursive algorithm begins (fc= 1) by setting 

ji(i) = i , ( 15b ) 


In order to choose jk+i { 0 , k = 1, . . . , rt — 1 , we consider as candidates the two stencils 

{ Xy t(0 _i, . . . , Zy t (0+k } or { Xy t(i) , . . . , Xy t (i)+fc+! } , (15c) 


which are obtained by adding a point to the left or right, respectively, of the previously 
determined stencil. We select the one in which tv is relatively smoother, i.t. the one in 


which the (/c + l)-th order divided difference is smaller in magnitude, thus 


V li\ = I “ 1 ’ if 

Jk+1 ^ ^ { jk{i ) , ot 


otherwise . 


'■ c jt(*)+ t ]l ^ I w [ x n(*)» • • • » x it(0+*+i] I * (15d) 
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Finally, we set j(i) = j n {i ) . 

For h sufficiently small, the algorithm (15) will always be able to determine an (n+1)- 
point stencil of smoothness between any two discontinuities. Therefore if each q n>i {: r; w) 
in (14b) is constructed with the aid of this adaptive-stencil type of interpolation, then 
wherever u>(x) is smooth, Q(x; w) will satisfy 

TV ( Q„(x ; u>) ) < TV (tu(x) ) + 0(h n+1 ) , (16) 


i.e. Q n (x ; w ) is essentially non-oscillatory. Henceforth, we will assume that Q n {x ; u;) is 
an ENO interpolating polynomial of degree n whose construction involves the adaptive- 
stencil algorithm (15). 

We now describe our reconstruction operator. We are initially given a discrete set of 
cell averages {©»,} , cell averages of a piecewise smooth function uj(x, y) , 


1 /•*■•+!/* rvi+i/2 . A . . 

©,-,- = - — — / / w{x,y)dydx 

Ax, A •'vy-i/a 


(17a) 


where Ax, Ay,- = a<, = (x,-+i/ 2 - £.- 1 / 2 ) (i/j+ 1/2 ~ l/y- 1 / 2 ) • F° r Vj+ 1/ 2 < V < ^-V 2 5 ^ e fi ne 
the primitive function W, (x) associated with w by 





1 rvj+1/2 

Ay, hj-i/. 


w{Z,y)dyd£ . 


(17b) 


Seeking a relationship between pointwise values of W,(x) and the discrete values {©,, } , 
we see immediately from the definitions (17a-b) that 


Ax,- uiij = Wj(x i+ i/i) - W,(x,-_i/ 2 ) , 


and we can therefore establish such a relationship at the cell interfaces. 


^■M=EAxt®ir ( 17c ) 

k—io 

Now, if we interpret the notation ©,(x) as the line average in y of u;(x,y) for a fixed 
x , then the definition (17b) clearly implies 


l V,+,/ 2 w i x >y) d y = ®i( x ) 

dx Ay,- V_ 1/a 


(17d) 


'Vy-1/3 


This suggests that if we approximate W, (x) by Q r (x; W,) , the ENO picewise polynomial 
of degree r which interpolates W, at the values given by (17c), we can then obtain an 
approximation of ©,(x) by defining the first step in our reconstruction procedure as 


R(x]iv) - — Q r {x]Wj) = 0,(x) . 


(I7e) 
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(17f) 


Then C ; (i) is a polynomial in x of degree r — 1 which satisfies 

Vj{x) = u )j{x) + 0(Ax r ) , 

wherever w, (x) is sufficiently smooth in x. 

If the procedure (17) is performed for all j , then we have a set of piecewise polynomials 
{^ (x)} , each of which is a high-order approximation to each t»,(x) . Clearly, from the 
definition (17d), the value of ©,(x) in x is equivalent in form to a one-dimensional cell 
average on the interval [yj+i/ 2 , yy- 1 / 2 ] ■ Therefore, for a fixed x, the remainder of our 
reconstruction procedure becomes equivalent to the one— dimensional method in [5], applied 
to the set {t? ; (x)} . 

For fixed x, we now treat the set {f>,(x)} as one-dimensional cell averages in y of a 
piecewise smooth function v(x, y ) , which we wish to reconstruct to high-order pointwise 
accuracy in y. Analogous to the method (17), we define another primitive function W (x, y) 
associated with v(x, y) by 

W (x, y) = f V v(x,y)dy , (18a) 

'Vo 

whose pointwise values we know at cell interfaces 

W{x,y j+1/2 ) = Ay* »*(*)• ( 18b ) 

k=jo 

Fitting the point values (18b) of W(y) with a piecewise polynomial Q r {y\W) of degree r 
by ENO interpolation, we can obtain a high-order pointwise approximation to v(x,y) in 
y by defining the second reconstruction step 

R(y ; f/(x)) = ^Q r (y,W) = p(x, y) , (18c) 

where, for fixed x , p(x, y) is a polynomial in y of degree r-1 that satisfies 

p(x,y) = t>(x,y) + 0(Ay r ) > (18d) 

wherever v(x,y) is sufficiently smooth. 

Noting the reconstruction definitions and the error relationships above, we can see that 
the values obtained in (18c) are the high-order pointwise approximations to w(x,y) which 
we desired from the initial cell averages 

iZ J (x,y; fi>) = R{y -,R{x;ib) ) = u>(x,y) + 0(h T ) , (19) 

We note here that the high-order relationship (19) cannot be achieved for r > 2 by simply 
applying the operator R to two u over lapped” one— dimensional stencils. 
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In addition to the high-order truncation error in regions of smoothness, we also note 
that each of the one-dimensional reconstruction operators (17e) and (18c) is essentially 
non-oscillatory, due to the nature of the interpolating polynomial Q r as given by (16). Fur- 
thermore, we note that R 2 is “conservative” , in the sense that the cell-averaging operator 
defined by the right-hand side of (17a) is its left-hand inverse, x.t. 


r*i+i/2 rvj+i/2 , =\ j j = 

/ / R {x,y\ w)dydx = Wij , 

» i /-a i 


A Xi A xjj Jyj- 1/3 


( 20 ) 


which is necessary in order that our numerical scheme (12) remain conservative. This 
property results directly from the various definitions in the reconstruction (17-18), and 
the fact that Q r is an interpolating polynomial. It is the adaptive-stencil algorithm (15) 
that enables this reconstruction (for sufficiently small h) to be high-order accurate on 
any domain where u>(x,y) is smooth, even if that region is near one in which w(x,y) 
is discontinuous. Furthermore, algorithm (15) is ultimately responsible for the adequate 
resolution of a discontinuity, near which the “jumps” in iZ J (x,y; ©) at cell boundaries 
become large relative to the mesh spacing. 

We further note that the error coefficient e(x,y) in (8), due to this reconstruction, 
becomes discontinuous at points where there is a change of orientation in the stencil of the 
associated interpolation. This discontinuity may occurr at critical points of the function 
and/or its derivatives. It is clear that when e(x, y) fails to be Lipschitz continuous at a 
point, the truncation error of the approximate flux difference in (10c) is only 0(h r+1 ) . We 
therefore expect the cumulative pointwise error due to N applications of the operator E^ 
to be only 0(h r ~ 1 ) at such points, but to remain 0(h r ) away from these points. Owing 
to the essentially non-oscillatory nature of E h , it is reasonable to expect the number of 
points at which e(x, y) fails to be Lipschitz continuous to remain bounded as h — ► 0. In 
this case, we see that the cumulative error of our numerical scheme is 0(h r ~ 1 ) in the Loo 
norm and 0(h r ) in the L\ norm. 


4. Systems of Conservation Laws 

In this section, we extend the scalar reconstruction procedure of Section 3 to solutions of 
hyperbolic systems of conservation laws. To this end, we now reconsider the IVP (1) , whose 
solution u is a vector of m components, as are the fluxes f(u) and g(u ) . We now wish to 
develop a vector reconstruction operator, denoted by R J , which will reconstruct a set {fijy} 
of vector— valued cell averages to high— order pointwise accuracy. Clearly, in this context, 
fLij denotes a vector whose components are the cell averages of the scalar components of 
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u. It would therefore seem natural, purely from an approximation theory viewpoint, to 
reconstruct the set by applying the scalar reconstruction R 2 in a component-wise 

fashion. However, this approach is valid only if we assume that u(x, y,f) remains free of 
discontinuities for the duration of a desired calculation. But in the more general case, 
when u(x,y, t) is a solution of a non-linearly coupled system of equations, such a solution 
can admit the collision of discontinuities of the same or of different families, as well as their 
collision with boundaries, t.g. solid walls. In the vicinity of such collisions, in the solution of 
more than one dependent variable, a component -wise reconstruction may develop spurious 
oscillations during this brief encounter which do not dissipate as the discontinuities then 
distance themselves from one another, if indeed they ever do. Numerical experiments to 
demonstrate this potential problem in one dimension can be found in [5], In the following, 
we describe an algorithm to reconstruct the vector-valued solution u from its cell averages 
which attempts to avoid this difficulty. 

We begin by considering the constant coefficient case of (1): /(u) = Au, y(u) = Bti , 
where A and B are constant m x m matrices; 

u t + Au, + Bu, = 0, (21a) 

u(x, y,0) = u°(x,y) . (21b) 

We also assume, for now, that our reconstruction takes place within the context of a 
Cartesian mesh. We note that the eigenvalues {A* (A)} of A and {A*(B)} of B are constant 
as are the right eigenvectors (r*(A)} , {r k {B)} and left eigenvectors {l k (A)} , {l k {B)} . We 
assume that these eigenvectors are suitably normalized so that 

/‘(A) • r>(A) = /*(B) -r’(B) = . (22a) 

If we define the Jfc th characteristic variables tv k (A) and w k (B) by the dot products 

to* (A) = l k {A ) u , u/*(B) = /*(£) u , k = 1 , 2 m , (22b) 

then it follows from (22a) that 

u = £ w»(A)r*(A) . = £ »*(B)r‘(B) . (22c) 

k=l *=1 

It is argued in [5] that, in the constant— coefficient case, it makes sense to use these 
characteristic variables (22b) in the reconstruction procedure, rather than u itself. This 
is due to the fact that, under the transformation (22b), the coupled system (2la) be- 
comes an uncoupled set of equations, thereby rendering any discontinuity in a particular 
characteristic variable “undetec table” by another. ( See [5] for details. ) 
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Therefore, in the case of a linear system (21a), we can describe our vector reconstruction 
as follows. Given {ft,,} , cell averages of a vector u, we begin by defining the cell averages 
of the characteristic variables in the x-direction by 

ft>‘(A) = l k {A) • Uij , for j fixed, all * , (23a) 

and then perform the scalar reconstruction given by (17e) on these averages. Using the 
result (22c), we can define the first step of our linear vector reconstruction procedure by 

R(x; ft) = £ R{x;W k {A))r k {A) = ft,(x) , (23b) 

jb=l 

the right-hand side of which is the vector-valued analogy of ft,-(x) in (17f). In analogy 
with the two-step procedure in Section 3, the reconstruction (23b) is performed for all j. 
For x fixed, we then proceed by approximating the “line-average” characteristic variables 
in the y-direction by the dot product 

w k (B) = l k {B) • ft, . (23c) 

The scalar reconstruction (18c) is applied to the values (23c) and, for a fixed x, we have a 
polynomial in y 

m 

R(y;ft(x)) = 53 R(y;© fc (S))r fc (B) , (23d) 

fc=i 

which completes our reconstruction for the linear case and we write 
R 2 (x,y; ft) = R( y ;R(x ; ft) ) 

= f; fi ( » ; l*(B) • f; *(*;/'(.!) -SKU) ) r*(B). (23e) 

*=i V p =i / 

We now wish to generalize the reconstruction procedure (23) to the case of a nonlinear 
system. In the nonlinear case of (21a), the matrices A(u) , B(u) are now functions of u, as 
are the eigenvalues {A*(A(u))} ,{A fc (B(u))} , and the eigenvectors {r*(A(u))} , {r*(R(u))} , 
{/*(A(u))} , and {/‘(B(u))} . Our extension will require the use of locally defined charac- 
teristic variables, in the following manner. In order to reconstruct u(x,y) on the region 
(x<-i/ 2 ,x, +1/2 ) x (y,-_i/ 2 ,y, •+!/*), we first derive a set of local average characteristic vari- 
ables {fl>*y(fti,)} , where n varies in the x-direction. We do this, for a fixed j, by computing 
dot products of / k (A(fi,,)) with the cell averages {ft„;} associated with intervals in the 
“immediate vicinity” of (x,_i/ 2 , x l+1 / 2 ) , 

©J,(fti,) = l k (A(Uij)) • fi„, , for n = * - q , . . . , * + q , (24a) 
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where q is the degree of the reconstruction polynomial. We then apply the scalar recon- 
struction operator R to this set {©*,•( tty)} of 2g+l variables in which the left eigenvector 
has been locally “frozen” at the I th location of the j th row of cells. This “local lineariza- 
tion” allows us to apply locally the linear vector reconstruction described in (23). The 
first step in our “nonlinear” vector reconstruction procedure then becomes 

R(x;tt) = £ R{x )r*(A(tt„)) = ©;(*) . (24b) 

Jb=l 

Upon performing (24b) for all j, then, for x fixed, we define a set of local “line— average” 
characteristic variables {^(ty)} in the y-direction : 

tD*(D,) = l k (B(Vj)) • 0 n , for n=j + q . (24c) 


We complete our “nonlinear” reconstruction by applying the scalar operator R to the set 
{©*(©,)} of local variables in (24c) which results in 


R(y ; ®(x)) 


£ 


Jb=l 




(24d) 


Thus, our vector reconstruction operator R 2 is a composition of (24b) and (24d) and, for 
a polynomial within a cell (t, j) , we write abstractly 


R J (x,y; fi)„- = E ^ y ; l*(B(«,(*)) •£ R(*;l'(A(i«)) .fi)r*(A(i<,)) r fc (B(e y (x)) , 

*=i \ p= i J 

(24e) 


and satisfies 

R 2 (x, y; ti(t)) = u(x,y, £) + 0(h r ) , (24f) 

wherever u is sufficiently smooth, and which we hope will avoid the oscillatory behavior 
associated with colliding discontinuities. 


5. CURVILINEAR CO-ORDINATES 

We now wish to generalize the spatial domain of solutions of the IVP (1). To this end 
let 

x = x(£,r?) , y = y{t,v) > ( 25a ) 

denote a smooth transformation from the physical x— y plane to the rectangular £— rj plane, 
with Jacobian determinant J given by 

J = x { y^ - y € x„ . (25b) 
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Let Cij denote a discrete region in physical space which is mapped into a rectangle 
(&- 1 / 2 , €«+ 1 / 2 ) x (r?y_i/ 2 ,r? J+1/2 ) by the transformation (25a). Then for any such region, our 
semi-discrete formulation for a weak solution of (la) remains identical in form to (2a) : 


S- = ft+l/2 - /»- 1/2, y (0 + 9 i,j+ 1/2(0 9 i,j- 1/2(0 

Ol Oij 

Under the transformation (25), we interpret the cell average fi,/(£) at time t as 
Qij(t) = — f 1+1/3 u(Z,r},t)J{t,T))dtdri , 

a ij *'£»- 1/3 

where a,y is the area of Cij • The fluxes f and g are given by 

/,.,(<) = P"" f («({,., /l.V.‘))dr, ; 

/*£«+ 1/3 ~ . 

ft , y+i/2(0 = / ’ 0 ) ’ 

■'Si- 1/3 

where 

F(u) = y n /(u) -x„y(u) , G'(u) = x € y(u) -y { /(u) , 


(26a) 

(26b) 

(26c) 

(26d) 

(26®) 


and /(u) and y(u) are the “known” flux vectors in (la). 

Having defined the necessary terms of our finite— volume formulation in a curvilinear 
co-ordinate system, we now wish to discretize our spatial and temporal operations. Again, 
with our Runge-Kutta time discretization, the extension is straight forward. As for the 
spatial operation, we note two basic modifications that must be made. 

The first alteration involves the reconstruction operator. Having clearly and consis- 
tently defined this procedure in terms of a rectangular mesh, we therefore would like to 
be able to perform the reconstruction procedure in £ — 27 space. Given a cell average fr, ; , 
as defined by (26b), we see that the quantity 6*, can be interpreted as a “cell average” 
of the function u(£,tj) J(£,»?) on the corresponding computational cell, having unit area. 
Therefore if we choose, instead, to reconstruct the set {a^- of cell averages which are 
“scaled” by their respective areas, we transform the reconstruction procedure to rectangu- 
lar co-ordinates, and can simply use the procedures outlined in Sections 3 and 4. In the 
case of vector reconstruction (Section 4), the various eigenvalues and eigenvectors required 
for this purpose are now the corresponding quantities of the matrices 


A = y n A - x n B , B = x i B-y i A. 


The polynomial p(£,»?) we construct by this approach will be a high-order approximation 
to the function u(£,r)) J{£, 27 ) , and therefore must be “re— scaled” by «/(£, 27 ) in order to 
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yield the pointwise values we need to approximate the fluxes (26c-d). Thus, we define our 
curvilinear reconstruction operator R 2 by 

fi) = R 2 {t>v;a*) , ( 27 ) 


with the vector analogy R 2 defined in an identical manner. 

As for our second modification, we note that when we desire accuracy that is higher 
than second— order, a piecewise linear interpretation of a two-dimensional grid will no 
longer suffice. This is due to the fact that once we require more than one point in the 
quadrature which will approximate the flux integrals (26c-d), we can no longer assume 
that the quantities {x^x n> y^y n } are constant along a given cell boundary. We must 
therefore assume that the mesh is “truly curvilinear” and account for any change in these 
grid metrics at each point required in the quadrature. It should also be noted that the 
curvature of cell boundaries also affects the lengths of cell interfaces, which become the 
curvilinear analogy of Ay, , Ax, in (12c-d). Our numerical flux can then be written 


7i+l/l,i(t) = £) Cfc F Rm [ R 2 {Zi+l/2 -0,Vk‘, %(*)) » R*{£i+l/2 + 0 > 

2 k= 1 




(28) 


where R 2 is given by (27), and |At?|<+t/>,j represents the arclength of the of the boundary of 
Cij along £ = &+ 1/2 • The description of g iJ+ i/ 2 (t) follows from symmetry. Now, assuming 
that the transformation (25a) is sufficiently smooth, the numerical scheme resulting from 
the curvilinear formulation (26) and the modifications mentioned above will be r —order 
accurate as defined by (4) in smooth regions and avoid oscillations near steep gradients. 

We make one further generalization. It so happens that there are a lot of applications of 
structured computational meshes for which a closed form transformation is not available. 
For ex am ple, a set of grid points may be initially generated as a solution of a system of 
differential equations, after which they may then be subjected to some smoothing operator, 
e.g. Laplacian. In such a case we do not have a set of equations (25a) from which to 
determine all the grid quantities necessary for the flux computation (28). 

However, given such a set of points, we might consider the equivalent of a locally de- 
fined set of transformation equations which are derived by polynomial approximation. By 
this we mean that each “grid line” through a set of points is approximated by polynomial 
interpolation, and all the necessary mesh quantities are calculated from these polyno- 
mials. As the grid metrics {x e , x„, y ( , y„} represent the components of outward normals 
at cell boundaries, a simple calculus argument will enable us to compute these quanti- 
ties at the desired quadrature points. Cell areas {cty} and the arclengths of cell faces 
{|A»7 1, •+!/*,,•, | A£|t ,^+ 1 / 2 } are also straight forward, once the four polynomials defining each 
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control volume are determined. An approximate Jacobian determinant, which we will de- 
note by J 9 , required for the curvilinear reconstruction procedure (27) is also obtainable 
by the following reasoning. If the transformation (25a) existed then clearly we would have 


the relationship 


dij — 


/■£•+ l/a /■■»/+ i/a 
l/a 


«J( £,»?)<*£<*? » 


and therefore we may interpret a local area as a rectangular “cell average” of J . Thus, 
having first determined approximations to the cell areas by polynomial interpolation of 
the grid points, we apply our scalar reconstruction algorithm to this set and define J' by 


•/'(£,*?) = R 2 {S,rn «) • 


We note here that, without a known transformation (25a) from which to mathematically 
determine “sufficient smoothness,” we can no longer make the claim concerning formal 
order of accuracy when this polynomial grid approximation is employed. We do, however, 
“test” this approach by implementing it in a numerical experiment rather than applying 
the exact grid quantities from the known transformation. The results of this test case are 
reported in Section 7. 


6. IMPLEMENTATION 

Before proceeding to numerical test cases, we find it prudent to make a few remarks 
concerning the implementation of some of the ideas in the preceding sections. As these 
notions are presented in a general setting, their implementation in a given problem is 
subject to individual interpretation and the suggestions made in this section are those 
which are employed by the author in the calculations presented in the next section. 

Our first topic concerns the first step in our scheme, namely the reconstruction pro- 
cedure. Though it is not necessary from an approximation theory viewpoint, we will 
implement the reconstruction in two separate stages, depending upon the orientation of 
the cell boundary along which we desire to approximate u. For the evaluation of the 
flux in (12c), the pointwise approximation of u at the Gauss points {y*} along 

x = x i+ i /2 is determined by 

u(x, +1 / 2 ±0,y*,t) « R( y* ; R{ x,- + i/j ± 0 ; fi(t) ) ) , (29a) 

and for the evaluation of h+i/i,j{t) ^ (l2d), we employ the composite operation (29a) in 
“reversed” form and achieve the required pointwise values by 

u(x t ,y,+x /2 ± 0,t) « R{x k ] R(y j+ i /2 ±0;fl(t)) ) . (29b) 
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Also pertaining to the reconstruction, we have taken the position that a high— order 
approximation must reflect available information. Therefore, near boundaries, we restrict 
our interpolation stencil to remain within the computational domain. It was found during 
numerical experimentation that problems could result from a high-order, one-sided inter- 
polation procedure. In particular, if a solution developed a shock which reflected from a 
solid wall, oscillatory behavior was noticed in the smooth region near the reflection point 
between the shock and the wall. Attempting to eliminate such numerical noise, one might 
suggest some sort of test for a “desirable reconstruction” which, when failed, will result in a 
local reduction in the order of interpolation. Owing to the recursive nature of the Newton 
interpolation procedure, such a test can be readily applied during the actual “building” of 
the polynomial. One such test simply checks for “over-shoots” or “under-shoots.” This 
is equivalent to requiring the endpoint values of the polynomial p,(x) which approximates 
u(x) on [xj_i/ 2 , £.+ 1 / 2 ] not to over-shoot a larger adjacent cell average, or to under-shoot 
a smaller one. We can so restrict p,(x) by requiring that it satisfy 


(&.+1 

- Ui )( Ui + 1 - Pi (x,+ 1/2)) > 0 

(30) 

(«, - 

«i-l) (Pi(*.--l/2) - tti-l) > 0. 

(31) 


We define any p,(x) which satisfies (30) to be a “desirable” reconstruction. If either 
of (30-31) are not satisfied, then the degree of p<(x) is reduced. Clearly, this restriction 
will never require any less than a locally linear reconstruction, as such a reconstruction is 
always locally monotone. 

As for the application of the boundary conditions themselves, such conditions as inflow 
or outflow are handled in a standard manner. As it turns out, a wall condition is relatively 
simple to handle. Because we actually approximate the solution on the boundary of each 
cell, those point values {u tt } along a solid wall are treated with an appropriate boundary 
condition and the numerical flux f w on the wall is simply 

h = f(WC(u w )) , (32) 

with WC{-) denoting that a wall condition has been applied. 


7. Numerical Experiments 

In the following, we present a few examples of numerically computed solutions using 
the scheme (12) as well as its extension to curvilinear co-ordinates. We have performed 
experiments on scalar equations in order to test for the computational order of accuracy, 
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and for our problems containing solid walls we solve the two-dimensional Euler equations 
of gas dynamics. 

In the tests involving the measurement of accuracy, we use r c to denote the “compu- 
tational order of accuracy.” This value is calculated by assuming a linear accumulation of 
error as given by (13). Also, in all of our results, when referring to a CFL coefficient, we 
imply its use in the “conventional” explicit time-step restriction, which in the scalar case 
is given by 

Ax Ay 

A ‘ = CFL max. ( |/'(«)| Ay + |^(»)| Ax) ’ 
and is easily extended to a hyperbolic system. 

EXAMPLE 1. In order to test our high-order ENO scheme (12) for its accuracy, we 
solve the two-dimensional linear advection equation 

u 4 + u r + u v = 0 , t > 0 , (33a) 

with initial data 

u(x, y,0) = i cos 7r (x + y) + | . (33b) 

The solution of (33) is 2-periodic in x and y for all time. By restricting our computational 
domain to -1 < x , y < 1 , we thereby make the boundaries 2-periodic also, effectively 
removing them from consideration. We note here that even though we are solving a 
linear equation, the high— order scheme (12) applied to (33a) is still non-linear, due to the 
adaptive stencil algorithm of the reconstruction operator R 2 . 

The exact solution of (33) can be easily calculated one— dimensionally in terms of the 
variable £ = x+y , and can be written 

u(x, y,t) = i cos 7r (x + y — 2t) + \ . (34) 

However, because we computationally solve (33) on a cartesian grid, our application is 
truly two-dimensional. We further emphasize the two— dimensionality of the numerical 
solution by discretizing the computational domain so that Ax ^ Ay . 

Since the solution (34) is smooth for all time, we apply our scheme for one period in 
time, ».«. t = 2.0, at a CFL of 2/3. The number of iterations required to do this on a 
given grid is high enough to expect a significant accumulation of error. We have performed 
this test on five consecutively refined meshes for the scheme (12) with orders of accuracy 
r = 1,2, 3,4. Though good results were obtained in all four cases, we are particularly 
interested in the “higher-order” cases r = 3 and r = 4 . These errors, along with then- 
corresponding computational orders of accuracy “r c ,” are presented in Table 1, and are 
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calculated with respect to the L and L\ norms. The exact and computational solutions 
were compared by cell-centered pointwise output. We see “one less order” of accuracy in 
the Loo norm for r = 4 . We therefore expect that there are points in the solution where the 
reconstruction stencil is discontinuous, as referred to in the closing paragraph of Section 
3. However, the third-order error for this particular problem is uniform. 

EXAMPLE 2. We now test a nonlinear equation, namely the two-dimensional Burgers 
equation 

+ + = °> * >0 ’ ( 35a ) 

again with initial data 

u(x, y,0) = icos7r(x + y) + i. (35b) 

We solve the IVP (35) on the same domain as the previous example, and again apply 
periodic boundary conditions. In this case, due to the non-linearity of equation (35a), 
gradients immediately begin to steepen upon the first time step, until a shock eventually 
forms at t = 2 / 7 T . We therefore apply the scheme using a CFL of 3/4 , up to t = 0.15 , 
when the solution remains smooth. Table 2 illustrates the accumulated errors for this 
test case for r = 3 and r = 4 . The exact solution is computed by using Newton-Raphson 
iterations to solve the characteristic relation 

u(x, y,t) = | cos 7r ( x + y — 2ut ) + | . 

In this case, we do notice the “drop” in the order of accuracy with respect to the L*, norm, 
in the third-order case. The solution was then computed to and past the point of shock 
formation, with no visible oscillatory behavior near the discontinuity. 

EXAMPLE 3. Our next numerical experiment involves the solution of the two- 
dimensional Euler equations of gas dynamics in a cartesian frame of reference, which we 
write in conservative form 


U t + F(U) X + G(U) y = 0 , 


(36a) 


where 



p 


pu 


pv 

u = 

pu 

pv 

, F(U) = 

pu 2 + P 

puv 

, G(£7) = 

puv 

pv 2 +P 


pE 


(pE + P) u 


( pE + P)v 


(36b) 
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Here p , P , and E are the density, pressure, and total specific energy, respectively, and u 
and v are the cartesian components of the velocity vector V . We close the system (36a) 
of four equations with the polytropic equation of state : 

P = ( 7-1 )p{E-\V 2 ), (36c) 

where 7 is the ratio of specific heats. 

We take our test case from [11]. A simple two-dimensional inlet with a step is installed 
onto a uniform cartesian grid, and the problem begins with a uniform Mach 3 flow directed 
towards the step, from left to right. As the length scales of the inlet’s configuration are 
identical to those in [11], we omit this description and and refer the interested reader to the 
cited reference. Scaling the time variable by the height of the inlet’s entrance, the solution 
reaches a steady state at approximately t = 12.0 . However, the structure of this steady 
state is relatively uninteresting, and therefore we compute the solution in a time accurate 
manner up to t = 4.0 , when the flow— field contains a complicated shock structure. 

At the problem’s outset, we assume that the inlet is filled with air which we model as 
an ideal gas, with 7 = 1 . 4 , with normalized initial free— stream conditions 

poo = 1*4 , Poo 1*0 , Uqo — 3.0 , Woo = 0 • (37) 

We implement the scheme (37) with a CFL of 0.8 on a 120 x 40 grid. For the flux 
computation in (12c-d), we approximate the required “Riemann fluxes” by the method 
developed by Roe [7], combined with an appropriate entropy correction (See e.g. [4].) 
The inflow boundary condition is specified by ( 12 ) and held fixed. Because the outflow is 
supersonic, the exit boundary condition has no effect on the flow, and therefore we assume 
all gradients to vanish at this boundary. At the walls, we apply the tangency condition 

V • h w = 0 , (38) 

where h w is the unit vector normal to the wall. The nature of this solution is such that 
the corner of the step is the center of a rarefaction fan, and hence is a singular point of the 
flow. We therefore apply a special treatment at this corner, as described in [11], in order 
to avoid large numerical errors generated in the neighborhood of this point which would 
hinder our qualitative comparison. 

Because second-order schemes are the current “state of the art,” we present second— 
and fourth-order accurate solutions at t = 4.0 in Figures la through Id, choosing density 
and Mach number as the variables of comparison. Both variables are plotted using thirty 
equally spaced contours. Sharper discontinuities are the most notable improvement in the 
fourth-order case, perhaps most notably the slip line emanating from the triple point near 
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the top wall. The Mach stem in this area is also more correct in its length and its position 
upstream. The weak shock from the comer of the step is also more pronounced in the 
fourth-order case, as is the other weaker slip line formed as this shock intersects the shock 
reflecting from the top of the step. This weaker slip line is virtually undetectable in the 
the second-order accurate solution on this computational mesh. 

EXAMPLE 4. Again, we solve the two-dimensional Euler equations of gas dynamics, 
this time as the solution pertains to the reflection of a moving shock wave from an inclined 
wall. The self-similar nature of such a solution lends itself to a rigorous analysis which is 
well documented in the literature ( e.g . [6]), and we therefore omit any general discussion 
of this phenomenon. These oblique shock reflections have been the subject of extensive 
experimental and computational research and the interested reader is referred to [1,2,12] 
and the references therein. 

Our problem begins with a plane shock, whose Mach number we denote by Ms , which is 
moving into still air towards a wall inclined by an angle 0 W to the direction of the shock s 
motion as shown in Figure 2a. The flow orientation is chosen to facilitate comparison 
with existing experimental interferograms. The problem becomes truly two-dimensional 
when the shock encounters the wall and forms a reflection whose structure can be quite 
complex. Analysis shows that the resulting similarity solution can be entirely determined 
by the parameters Ms and 0 W . 

We examine two cases, both of which are double Mach reflections. The wall angle 0 W 
is 40 degrees in both cases, and the shock Mach numbers we examine are 2.87 and 3.72. 
This type of reflection exhibits a complex structure containing shock diffractions and slip 
lines, and is particularly demanding of any computational algorithm. For this reason, these 
computations are most commonly performed with the use of a self-similar transformation 
which effectively removes the time dependence of the solution. The resulting equations 
that are then solved resemble the steady Euler equations with the addition of source terms. 
However, because we wish to examine the temporal as well as the spatial accuracy of our 
scheme, we choose to compute our solutions in a time-accurate manner. 

In addition to the demanding nature of the solution itself, such computations are also 
made difficult by geometric concerns, largely due to the presence of a sharp corner. Though 
there is no way to rid ourselves of the corner itself, we attempt to mitigate its presence 
by using a curvilinear grid transformation. A portion of our particular mesh is shown 
in Figure 2b, and is generated by a Schwarz-Chrisoffel transformation. We could use 
this transformation to derive all of the necessary grid quantities referred to in Section 5. 
However, we would like to test our scheme in its most general form. Therefore, given 
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a collection of points generated by this transformation, we calculate our all of our mesh 
variables from the approximate grid lines we generate by polynomial interpolation. 

The equations we solve are the time— dependent Euler equations of gas dynamics in 
general curvilinear co-ordinates, which we write 

U t + F(U) e + G{U) n = 0 , (39a) 

where 

F = y n F-x r ,G, G = x i G-y i F , (39b) 

and F and G are the cartesian flux vectors in (36b). Our initial conditions are the two 
constant states U 0 and U\ which determine the desired plane shock, oriented as in Figure 
2a. We normalize these conditions with respect to the still-air initial state Uq . These 
conditions then are 

p 0 - 1.0 , Po = 1.0 , «o = Vo = 0 . (40) 

Using these conditions and the shock Mach number Ms , we can then determine the state 
Ui by means of the Rankine-Hugoniot jump condition. These initial states axe then con- 
servatively interpolated onto the computational mesh. 

Provided we ensure that the entire reflection remains within the limits of the grid, the 
boundary conditions are relatively simple. The initial conditions Uq ,U\ are applied on the 
left and right boundaries, while the only concern at the fax-field is adequately accounting 
for the movement of the plane shock. On the wall, the tangency condition (38) is imposed. 
We apply our scheme in its curvilinear version to the two test cases for 400 time steps 
using a CFL of 0.8 on a 180 x 40 grid. Our purpose is not only to compare our numerical 
solutions by their order of accuracy, but also to compare each of them with experimental 
results, which we obtain from [3]. 

Figure 3a is an inteferogram resulting from an experiment designed to photographically 
exhibit the density structure for the case Ms = 2.87 . (The alphabetical labeling of this 
picture is not relevant to our presentation.) Density contour plots for the second- and 
fourth-order numerical solutions axe compared in Figures 3b and 3c. Overall, crisper 
discontinuities are observed in the fourth-order solution. The shock structure itself is also 
more correct in Figure 3c, in that the more perpendicular orientation of the incident Mach 
stem with respect to the wall is more in line with the experimental observation. The 
“toeing out” of this Mach stem in Figure 3b appears to be due the poorer resolution of 
the contact discontinuity emanating from the primary triple point. 

In Figure 3d, we plot the fourth-order density solution on the wall in order to make a 
comparison with experimental measurements. The x — axis is scaled by the distance L from 
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the incident Mach stem on the wall to the corner, the Mach stem location being x — 0, 
and the the comer at x = 1 . Overall agreement with the experimental data is as good if 
not better than the similar comparison in [3] in which the numerical solution was achieved 
in a self-similar fashion and on a much finer grid. 

In Figures 4a through 4d, we have the analogous results for the case M s = 3.72 . In 
Figure 4b, the resolution of the contact discontinuity is much worse than in Figure 3b, 
causing a serious error in the formation of the incident Mach stem. There is no such 
problem, however, in the fourth-order case (Figure 4c) where the resolution of the contact 
discontinuity is excellent. Finally, the density wall plot in Figure 4d appears to be an 
improvement over the similar result in [3]. 


8. Concluding Remarks 

The usefulness of shock-capturing schemes which are of formal high-order accuracy 
in more than one spatial dimension appears certain. For instance, the results in EXAM- 
PLE 3 suggest that these types of schemes might provide better results than are currently 
available in scientific areas where the resolution of weak waves is crucial, t.g. acoustics. 
EXAMPLE 4 illustrates the effect higher-order accuracy can have in the more accurate 
resolution of unsteady compressible flow solutions which develop complicated shock struc- 
tures. Furthermore, when compressible flow solutions are also required to be viscous, 
these schemes can play a major role in the computation of such flows, in which there are 
regions containing many smooth local extrema. In addition, when such flows also develop 
shocks, c.ff. shock-turbulence interactions, the need for these newly-developed schemes is 
undeniable. 
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TABLE 1 

Solution Error for IVP (8) 
t = 2.0 CFL= | 


Loo error error 


Grid 

r = 3 

r c 

r = 4 

r c 

r — 3 

r c 

r = 4 

r c 

8 x 12 

1.559 E-l 


8.764 E-2 


3.890 E-l 


2.078 E-l 


16 x 24 

2.331 E-2 

2.74 

1.021 E-2 

3.10 

5.678 E-2 

2.78 

1.605 E-2 

3.69 

32 x 48 

3.002 E-3 

2.97 

1.239 E-3 

3.04 

7.627 E-3 

2.97 

1.764 E-3 

3.20 

64 x 96 

3.749 E-4 

3.00 

1.390 E-4 

3.16 

9.130 E-4 

2.99 

1.237 E-4 

3.82 

128 x 192 

4.666 E-5 

3.01 

1.544 E-5 

3.17 

1.142 E-4 

3.00 

8.658 E-6 

3.84 


TABLE 2 

Solution Error for IVP (10) 
t = 0.15 CFL= f 

Loo error Li error 


Grid 

r = 3 

r c 

r = 4 

r e 

r = 3 

r e 

r — 4 

r e 

8 x 12 

3.765 E-2 


2.632 E-2 


4.938 E-2 


3.010 E-2 


16 x 24 

9.549 E-3 

1.98 

4.373 E-3 

2.59 

8.661 E-3 

2.51 

3.834 E-3 

2.97 

32 x 48 

2.111 E-3 

2.18 

4.192 E-4 

3.38 

1.237 E-3 

2.81 

3.245 E-4 

3.56 

64 X 96 

4.264 E-4 

2.31 

4.374 E-4 

3.26 

1.844 E-4 

2.75 

2.497 E-5 

3.70 

128 x 192 

8.988 E-5 

2.25 

4.361 E-6 

3.33 

2.783 E-5 

2.73 

1.968 E-6 

3.67 
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